#ifndef SPECIES1D_H
#define SPECIES1D_H

#include <vector>
#include <string>

#include "../Species/Species.h"

using namespace std;

class Species1D : public Species
{
public:
    //constructor
    Species1D(PicParams* params_, ParallelMPI* pmpi_, int ispec);

    //destructor
    ~Species1D();

    //calculate the particle dynamics (interpolation, pusher, projection)
    //only caculate the number density, no electric currents
    virtual void move(PhysicalField* fields);

    //project particles to grid points, to calculate charge density
    virtual void project(PhysicalField* fields);

    //charged particle: interpolate fields at particle position, and move particle
    void push(CellParticles<1>* cell_particles, const LocalFields E_cell_fields[2], const LocalFields B_cell_fields[2]);

    //neutral particle: interpolate fields at particle position, and move particle
    void push_neutral(CellParticles<1>* cell_particles, const LocalFields E_cell_fields[2], const LocalFields B_cell_fields[2]);

    //decide if particles hit boundary, cross MPI boundary, to move to another cell, or still stay in the cell
    virtual void sort_particles();

    //After MPI exchange particles, delete useless particles
    virtual void correct_particles();

    inline void interpolate(const LocalFields cell_fields[2], LocalFields& local_field, const double position[1])
    {
        //Normalized particle position
        double xpn = position[0] * params->cell_length_inv[0];

        //Indexes of the central nodes
        int ip_ = floor(xpn);

        //Declaration and calculation of the coefficient for interpolation
        double delta;
        double coeffxp[2];

        delta   = xpn - (double)ip_;
        coeffxp[0] = 1.0 - delta;
        coeffxp[1] = delta;

        local_field.x = 0.0;
        local_field.y = 0.0;
        local_field.z = 0.0;

        for (int iloc=0 ; iloc<2 ; iloc++)
        {
            local_field.x += coeffxp[iloc] * cell_fields[iloc].x;
            local_field.y += coeffxp[iloc] * cell_fields[iloc].y;
            local_field.z += coeffxp[iloc] * cell_fields[iloc].z;
        }
    };

    virtual int get_number_of_particles()
    {
        int n = 0;
        for(int i = 0; i < particles_.size(); i++)
        {
            n += particles_[i]->size();
        }
        return n;
    }

    void init_particles();

    vector< CellParticles<1>* > particles_;

    vector< Particle<1> > psi_particles_;

    vector< Particle<1> > exchange_particles_;

    PartBoundCond<1> part_bound_cond;


};



#endif
